# get data
load(file="scripts/cERGM_data.RData")

library(statnet)
library(MLmetrics)
library(xtable)
library(tc.ergmterms)

# there are cases in 1962 marked as 1963 and in 1971 marked as 1972
scc1[3364,11] <- 1962
scc1[4648,11] <- 1971


year.total <- scc1[,11]-1936 #-1936 => the first year 1937 is 1, 1938 is 2 aso, term column 

set.seed(12123)


# beginning of big loop. takes a few weeks to run. suggest to split up the loop and let it run parallel.
for(t in 14:79){  # 14 corresponds to 1950 and 79 to 2015
print(t+1936)
tt<- t+1936
rowlasttt <- max(which(scc1[,11] == tt))
idt <- scc1[rowlasttt, 55]  


cases <- max(which(year.total==t)) # 14 corresponds to 1950

# simple assignment of time periods to cases
years <- year.total[1:cases] 

# going to need a sender time matrix covariate
sender.time <- matrix(years,length(years),length(years),byrow=F)
year <- matrix(years, length(years),length(years),byrow=F)

# extract the network up to time t
AM <- adjacency.matrix[which(years <= t),which(years <= t)]

net.t <- network(AM) 

#calculate outdegree
o.degree <- rowSums(AM)

# fix the outdegrees of time t as 0
last.t <- which(years==t)

o.degree[last.t] <- 0

# indicate outdegree and term as a node attribute
net.t <- set.vertex.attribute(net.t,"o.degree", o.degree)
net.t <- set.vertex.attribute(net.t,"term", years)


# subset MQ matrix
mq.t <- mq.matrix[which(years <= t),which(years <= t)]
# subset same issue matrix
same.issue.area.t <- same.issue.area[which(years <= t),which(years <= t)]
# subset year diff matrix
year.diff.t <- year.diff.matrix[which(years <= t),which(years <= t)]
# subset year diff matrix square
year.diff.square.t <- year.diff.t^2
# subset sender time matrix
sender.time.t <- sender.time[which(years <= t),which(years <= t)]
# subset year matrix
year.t <- year[which(years <= t),which(years <= t)]
# subset same opinion writer matrix
same.opinion.writer.t <- same.opinion.writer[which(years <= t),which(years <= t)]

#### set vertex attributes
# same issue area
net.t <- set.vertex.attribute(net.t,"SameIssueArea", scc1[which(years <= t),41])
# abs diff of MQ score
net.t <- set.vertex.attribute(net.t,"AbsDiffMQscores", scc1[which(years <= t),65])
# number justices that voted for the case
net.t <- set.vertex.attribute(net.t,"NumberJusticesPro", scc1[which(years <= t),52])
# overruled covariate
net.t <- set.vertex.attribute(net.t,"Overruled", Overruled.matrix[which(years <= t),idt])
# sender time
net.t <- set.vertex.attribute(net.t,"sender.time", sender.time.t[,1])
# Majority Opinion Writer
net.t <- set.vertex.attribute(net.t,"MajOpWriter", scc1[which(years <= t),49])

not.fixed <- network(1*(sender.time == t))

# get statistics
obs <- summary(net.t~edges+ mutual + nodeicov("o.degree")+ difftransties("term")+gwidegree(1, fixed=TRUE)+ 
                 dgwesp(0.15, fixed=TRUE, type="OSP") +
                 edgecov(mq.t)+
                 edgecov(year.diff.t)+ edgecov(year.diff.square.t )+ nodeicov('AbsDiffMQscores')+ nodeicov('NumberJusticesPro')+ 
                 nodeicov('Overruled')+edgecov(same.opinion.writer.t) + nodeifactor('MajOpWriter')+ edgecov(same.issue.area.t)  )
obs


san.net <- san(net.t~edges+ mutual + nodeicov("o.degree")+ difftransties("term")+ 
                 gwidegree(1, fixed=TRUE)+ 
                 dgwesp(0.15, fixed=TRUE, type="OSP") +
                 edgecov(mq.t)+
                 edgecov(year.diff.t)+ edgecov(year.diff.square.t )+ nodeicov('AbsDiffMQscores')+ nodeicov('NumberJusticesPro')+ 
                 nodeicov('Overruled')+edgecov(same.opinion.writer.t) + nodeifactor('MajOpWriter')+ 
                 edgecov(same.issue.area.t), constraints=~fixallbut(not.fixed) , target.stats = obs)


## stats of san network
summary(san.net~edges+ mutual + nodeicov("o.degree")+ difftransties("term")+ 
          gwidegree(1, fixed=TRUE)+ 
          dgwesp(0.15, fixed=TRUE, type="OSP") +
          edgecov(mq.t)+
          edgecov(year.diff.t)+ edgecov(year.diff.square.t )+ nodeicov('AbsDiffMQscores')+ nodeicov('NumberJusticesPro')+ 
          nodeicov('Overruled')+edgecov(same.opinion.writer.t) + nodeifactor('MajOpWriter')+ 
          edgecov(same.issue.area.t)  )


# get MPLE of SAN network

mple.san <- ergm(san.net~edges+ mutual + nodeicov("o.degree")+ difftransties("term")+
                   gwidegree(1, fixed=TRUE)+ dgwesp(0.15, fixed=TRUE, type="OSP") +
                   edgecov(mq.t)+
                   edgecov(year.diff.t)+ edgecov(year.diff.square.t )+ 
                   nodeicov('AbsDiffMQscores')+ nodeicov('NumberJusticesPro')+ 
                   nodeicov('Overruled')+
                   edgecov(same.opinion.writer.t)  + edgecov(same.issue.area.t),
                 constraints=~fixallbut(not.fixed), estimate="MPLE"  ) #
li <- coef(mple.san)



# fitting glm
model <- ergm(net.t~edges+ mutual + nodeicov("o.degree")+ difftransties("term")+gwidegree(1, fixed=TRUE)+ dgwesp(0.15, fixed=TRUE, type="OSP") +
                edgecov(mq.t)+
                edgecov(year.diff.t)+ edgecov(year.diff.square.t )+ 
                nodeicov('AbsDiffMQscores')+ nodeicov('NumberJusticesPro')+ 
                nodeicov('Overruled')+
                edgecov(same.opinion.writer.t)  + edgecov(same.issue.area.t),
              constraints=~fixallbut(not.fixed), control=control.ergm(MCMC.samplesize=15000,MCMLE.maxit=20, init=li)) #
summary(model)
aicm <- AIC(model)
bicm <- BIC(model)

#############################################
# Independence model
#############################################
ind <- ergm(net.t~edges+ 
              edgecov(mq.t)+
              edgecov(year.diff.t)+ edgecov(year.diff.square.t )+ 
              nodeicov('AbsDiffMQscores')+ nodeicov('NumberJusticesPro')+ 
              nodeicov('Overruled')+
              edgecov(same.opinion.writer.t)  + edgecov(same.issue.area.t),
              constraints=~fixallbut(not.fixed), estimate="MLE") #
summary(ind)   
aici <- AIC(ind)
bici <- BIC(ind)

li <- list(coef(model),  coef(summary(model))[,5], coef(summary(model))[,2], aicm, bicm,
          coef(ind),  coef(summary(ind))[,5], coef(summary(ind))[,2], aici, bici ) 
names(li) <- c("Coefficients", "P-values", "StdErrors", "AICm", "BICm",
              "CoefficientsI", "P-valuesI", "StdErrorsI", "AICi", "BICi")

yy <- tt-1900
foo <- paste('li', yy, sep ='')

assign(foo, li)

save(list=foo, file=paste(tt, "san.RData", sep="_"))
}

#################################################
## load all R-spaces that were created in the loop above

load(file="1950_san.RData")
load(file="1951_san.RData")
load(file="1952_san.RData")
load(file="1953_san.RData")
load(file="1954_san.RData")
load(file="1955_san.RData")
load(file="1956_san.RData")
load(file="1957_san.RData")
load(file="1958_san.RData")
load(file="1959_san.RData")
load(file="1960_san.RData")
load(file="1961_san.RData")
load(file="1962_san.RData")
load(file="1963_san.RData")
load(file="1964_san.RData")
load(file="1965_san.RData")
load(file="1966_san.RData")
load(file="1967_san.RData")
load(file="1968_san.RData")
load(file="1969_san.RData")
load(file="1970_san.RData")
load(file="1971_san.RData")
load(file="1972_san.RData")
load(file="1973_san.RData")
load(file="1974_san.RData")
load(file="1975_san.RData")
load(file="1976_san.RData")
load(file="1977_san.RData")
load(file="1978_san.RData")
load(file="1979_san.RData")
load(file="1980_san.RData")
load(file="1981_san.RData")
load(file="1982_san.RData")
load(file="1983_san.RData")
load(file="1984_san.RData")
load(file="1985_san.RData")
load(file="1986_san.RData")
load(file="1987_san.RData")
load(file="1988_san.RData")
load(file="1989_san.RData")
load(file="1990_san.RData")
load(file="1991_san.RData")
load(file="1992_san.RData")
load(file="1993_san.RData")
load(file="1994_san.RData")
load(file="1995_san.RData")
load(file="1996_san.RData")
load(file="1997_san.RData")
load(file="1998_san.RData")
load(file="1999_san.RData")
load(file="2000_san.RData")
load(file="2001_san.RData")
load(file="2002_san.RData")
load(file="2003_san.RData")
load(file="2004_san.RData")
load(file="2005_san.RData")
load(file="2006_san.RData")
load(file="2007_san.RData")
load(file="2008_san.RData")
load(file="2009_san.RData")
load(file="2010_san.RData")
load(file="2011_san.RData")
load(file="2012_san.RData")
load(file="2013_san.RData")
load(file="2014_san.RData")
load(file="2015_san.RData")


##########################################
## save an R-space with all loaded R-spaces.

save.image(file="Model_fit_results.RData")
